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Abstract 

In this paper, we present an analytical solution to nonlocal continuum electrostatics for an arbitrary 
charge distribution in a spherical solute. Our approach relies on two key steps: (1) re-formulating the 
PDE problem using boundary-integral equations, and (2) diagonalizing the boundary-integral operators 
using the fact their eigenfunctions are the surface spherical harmonics. To introduce this uncommon 
approach for analytical calculations in separable geometries, we rederive Kirkwood's classic results for 
a protein surrounded concentrically by a pure-water ion-exclusion layer and then a dilute electrolyte 
(modeled with the linearized Poisson-Boltzmann equation). Our main result, however, is an analytical 
method for calculating the reaction potential in a protein embedded in a nonlocal-dielectric solvent, 
the Lorentz model studied by Dogonadze and Kornyshev. The analytical method enables biophysicists 
to study the new nonlocal theory in a simple, computationally fast way; an open-source MATLAB 
implementation is included as supplemental information. 

1 Introduction 

One of the long-standing challenges in molecular biophysics is the development of accurate, yet simple models 
for the influence of biological fluids (aqueous solutions composed of water and dissolved ions) on biological 
molecules such as proteins and DNA. Atomistic simulations that include explicit water molecules, such 
as molecular dynamics (MD), provide the most detailed molecular understanding that is widely accessible 
without specialized computational resources. However, these simulations come at two prices: first, MD 
simulations can require many hundreds of compute hours, most of which are spent on the thousands of water 
molecules whose individual behaviors are not of primary relevance; second, practitioners must understand 
numerous subtleties about simulation protocols and the parameters associated with the physical models 
(force fields). Implicit-solvent models replace the explicit water molecules with an approximation to the 
theoretically rigorous potential of mean force (PMF) [U, creating the possibility of simulating molecular 
behavior accurately but orders of magnitude faster, and with fewer statistical uncertainties. Unfortunately, 
the statistical mechanical derivation of the PMF is not constructive, in the sense that the derivation does 
not provide a general PMF suitable for all molecular solutes. Instead, one must guess a functional form, 
such as the Poisson equation for the electrostatic interactions between solvent and solute, find the optimal 
parameters, and then test its fit against real data (both experiment and more accurate theories such as MD). 

Of course, evaluation of an implicit-solvent model is greatly accelerated if it can be solved easily and 
rapidly on relevant, non-trivial problems. With the advent of fast computers, one reasonable option is to 
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make numerical software implementing the new model freely available online [HI [37] ■ Another option is to 
provide analytical solutions for tractable geometries. Spheres are frequently used for continuum electrostatic 
modeling, because exact results can be obtained using spherical harmonics and the method of separation 
of variables [351 132] ■ Kirkwood's classic solution for a spherical protein embedded in a dilute electrolyte 
represents the best-known example [35], and demonstrates this conceptually simple approach. One merely 
writes down spherical-harmonic expansions and matches expansion coefficients using the known boundary 
conditions. Even though proteins obviously have complicated shapes, analysis of spherical geometries can 
offer insights into problems such as pA' a predictions |26j . redox potentials |59j . strategies for optimizing 
molecular binding [33] , and fast analytical models such as Generalized Born [5T| . 

However, Kirkwood's work also demonstrates a difficulty with the approach: as one adds detail to 
the model — in Kirkwood's case, an ion-exclusion layer outside the protein — calculations become onerously 
complex very quickly. Every additional layer or unknown function introduces another set of expansions that 
need to be matched, and manual algebraic manipulation for the desired expansion coefficients essentially 
entails solving a linear system of equations, so that the number of operations grows cubically with the number 
of equations. In addition, modeling the linearized Poissom-Boltzmann equation in the solvent necessitated 
the introduction of a set of polynomials for the radial coordinate because the standard Bessel functions were 
unsuitable [35]; more than sixty years passed before the relationship between Kirkwood's polynomials and 
the Bessel functions was established, allowing at the end a substantial simplification [42] . 

In this paper, we present an alternative strategy for obtaining analytical solutions in separable geometries. 
The first step is to transform the given system of partial-differential equations (PDEs) into one of boundary- 
integral equations (BIEs) [B], so that the unknowns are no longer functions defined over three-dimensional 
regions of space, but instead functions defined on two-dimensional boundaries. Second, the boundary-integral 
operators are diagonalized using the appropriate harmonics |181 136] . This allows a mode-by-mode calculation 
of the unknown functions on the boundary in terms of the appropriate surface harmonics — in contrast 
to matched-expansion approaches that employ solid harmonics. To demonstrate the BIE-eigenfunction 
approach, we solve the Kirkwood problem (a spherical protein embedded in a dilute electrolyte, with a thin 
ion-exclusion or Stern layer [35] ) and derive the full solution to the more recent nonlocal-diclcctric model of 
Dogonadze and Kornyshev [21] [38] [54] . 

The nonlocal model was originally developed to address one of the key shortcomings of macroscopic 
continuum theories for molecular solvation: the fact that the solvent molecules (usually water) are not 
infinitesimally small compared to length scales of interest, e.g., small ions [T5] [53] and proteins [35] ■ Un- 
fortunately, nonlocal response means that even the simplest form of the nonlocal model, called the Lorentz 
nonlocal theory [7J, leads to an integrodifferential Poisson equation, which is difficult to solve analytically 
or even numerically. Consequently, to date the only analytically solved geometries for the Lorentz nonlocal 
model have been the sphere with central charge [16j [56] and the charge near a half-space [46j [47j [45] , and 
no numerical algorithms for the original nonlocal model in arbitrary geometries were ever presented. 

Very recently, however, Hildebrandt and collaborators derived several mathematical reformulations to 
render the Lorentz nonlocal electrostatic model tractable both analytically and computationally [5^] HH1 
l30l [57] . The first major step was reformulating the nonlocal integrodifferential Poisson problem in one 
unknown variable, the electrostatic potential 9?(r), as a pair of coupled, purely local PDEs with two unknown 
variables throughout space (</?(r) and an additional auxiliary potential) |29j . Similar reformulations of 
nonlocal continuum theory were obtained independently in other areas of physics [33J [22] . Following this 
reformulation, Green's theorem and double reciprocity can be used to transform the coupled PDE system 
into a purely boundary-integral-equation (BIE) representation of the nonlocal model [30] [23] . 

In principle, both the local-formulation PDE problem and the purely BIE method are solved problems 
numerically, in the sense that asymptotically optimal (linear-scaling) numerical algorithms exist [9] [17] [41] 
[5j[52[TT]. However, even "fast solvers" can require an hour or more of computation, and therefore analytical 
solutions of non-trivial problems still hold significant value in this relatively early stage of testing nonlocal 
electrostatics of molecular solvation. One application of analytical methods is to obtain qualitative insight 
into the differences between nonlocal and local models using visualization: analytical methods allow rapid 
calculations of the reaction potential induced throughout a model geometry by a chemical group in the 
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protein, e.g. an amino acid side chain. Another application of analytical methods is to obtain quantitative 
information that may help to determine model parameters. For example, the nonlocal model includes 
an additional parameter beyond those of the standard local model. This parameter, denoted by A, is an 
effective length scale that captures water's transition from behaving like a low-dielectric material at short 
length scales to more familiar high-dielectric, bulk-like behavior at longer length scales. Parameterization 
requires extensive simulation and testing, and fast calculations aid significantly. 

To support the development and testing of nonlocal electrostatic models for biomolecule solvation, we 
present here the nonlocal-model analogue of Kirkwood's result: namely, an analytical approach for the elec- 
trostatic solvation free energy of an arbitrary charge distribution in a spherical solute embedded in a solvent 
modeled as a Lorentz nonlocal dielectric. Kirkwood's classic work continues to have impact decades after 
the advent of numerical simulations of the continuum electrostatic model |55[ 1261 151) , and the present work 
significantly enlarges the scope of nonlocal problems that can be studied analytically. We note that mobile 
ions such as sodium and potassium play crucial physiological roles and that the present work addresses only 
pure water solvent. However, the nonlocal theory can be extended easily to linearized Poisson-Boltzmann 
treatment of physiological electrolyte solutions [33], and these extensions are the subject of ongoing work. 

The remainder of the paper is organized as follows: the next section describes the local and nonlocal 
models, their reformulation as systems of boundary-integral equations, and the eigendecompositions of the 
associated boundary-integral operators. In Section [3] we introduce our BIE-cigenfunction strategy by red- 
eriving the solution to Kirkwood's problem, and then apply the strategy to solve the nonlocal problem. In 
Section 21 we present several applications of the analytical solution, which illuminate important differences 
between local and nonlocal electrostatics, including the choice of solute dielectric constant and the sensitivity 
of the nonlocal results to the solvent length-scale parameter A. The paper concludes in Section [5] with a brief 
summary and discussion. 

2 Background 



^water 




Figure 1: Diagram of the two continuum electrostatic models to be solved analytically, (a) Kirkwood's 
problem |35| . (b) Nonlocal- response model in a pure- water solvent. 
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2.1 Kirkwood's Local- Response Electrostatic Model 

Figure [TJa) is an illustration of the local-response model under consideration. We assume that the solute 
region I is a sphere of radius b, which is centered at the origin, and that the solute is at infinite dilution in 
a dilute aqueous electrolyte solvent. The solute charge distribution p(r) is modeled as a set of Q discrete 
point charges contained within the sphere, the ith of which has value qi and is situated at (r,;, 9i, </>j). The 
solute is treated as a homogeneous local-response dielectric with relative permittivity e pro tcin, i-e. inside the 
protein, the constitutive relation between the displacement and electric field is 

Di(r) = e pr otemeoEi(r) (1) 

where as usual E(r) = — V<y3(r) with ip the electrostatic potential. Substituting this constitutive relation 
into Gauss's law for dielectrics 

V • Di(r) = p(p), (2) 
we see the electrostatic potential in region 1 satisfies the familiar Poisson equation 

VV« = -A. (3 ) 

cin 

In a thin solvent layer surrounding the protein, we have water but no mobile ions; assuming that they are 
point charges in hard spheres of radius d, the ion density must be zero for ||r|| < b + d. Consequently, in this 
region (labeled II in Figure Ha)) the potential satisfies a Laplace equation and we assume the permittivity 
is just that of pure water e wa ter ~ 80. Standard boundary conditions hold at the protein-solvent interface 
defined by ||r|| = b, namely the continuity of the potential and the normal component of the displacement 
field: 

<pi(*b) = <m( T t) (4) 

n-D I (r b -)=n-D„(r+). (5) 
For local-response dielectrics, Eq. [5] reduces to the familiar 

c protem — c watcr • V u / 

where the superscripts — and + denote the interior (solute) and exterior (solvent) regions, respectively, and 
the normal direction n points outward from region I to region II. 

Outside this ion-exclusion layer, the mobile ions are assumed to redistribute such that at any point r, 
the net charge density is the sum of the Boltzmann- weighted ion densities (i.e., neglecting the ion sizes and 
correlations between them). This leads to the nonlinear Poisson-Boltzmann equation, which here we simplify 
by linearization, i.e. the potential in region III satisfies the linearized Poisson-Boltzmann equation (LPBE) 

VVm(r) = K 2 vra(r) (7) 

where k is the inverse Debye screening length; for physiological solutions, n « 8 A. The electrolyte is also 
assumed to have relative permittivity e wa t er , and so the boundary conditions at the ion-exclusion boundary 
||r|| = a are 

^ii(r-) = wn(r+) (8) 



dn dn 



(9) 



Kirkwood solved the above problem for the potential using matched expansions in the solid spherical 
harmonics [35] . Here, we show that an alternative is to use the surface harmonics for the BIE formulation of 
this problem, which may be derived as follows. For a point r in one of these regions, Green's representation 
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theorem allows the potential at r to be written in terms of the potential and its normal derivative at the 
surface or surfaces the bound the region [32] [33j [58| [5] . In region I, for example, 



(r) 



dG L (r,r') 
dn 



Vl (v')dA' - [ G L (r, r ')^ldA' + f G L (r, v')p(v')dV, 

Jb l ~' n J region I 



(10) 



where the subscript b denotes the spherical boundary ||r|| =b, G L (r, r') 



is the free-space Green's 



4?r| r— r'| 

function for the Laplace equation, and the third term on the right-hand side represents the Coulomb potential 
induced by the solute charge distribution. Writing similar expressions for the potential in regions II and III, 
and taking careful limits as the field points approach these bounding surfaces, we obtain a system of four 
boundary-integral equations for the four unknown functions (the potential and normal derivative on the two 
boundaries). The complete derivation is presented elsewhere [5], but the final system may be written as 
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(11) 



Here, we have introduced a short-hand operator notation in which / denotes the identity operator, V denotes 
a single- layer potential operator (the second term on the right-hand side of Eq. [TU)) and K denotes a double- 
layer potential (the first term on the right-hand side of Eq. m?]) ; the superscripts L and Y denote the 
Laplace or linearized Poisson-Boltzmann (Yukawa) Green's function; and the subscript pair 6, a denotes the 
"source" surface (a) and the "destination" surface (b). The identity-operator terms arise from singularities 
in the double-layer potential. 



2.2 Nonlocal-Response Electrostatic Model 

Figure (Hb) is an illustration of the nonlocal-response model. As in the local-response problem, we assume 
a spherical solute of radius b, centered at the origin, with Q discrete point charges as the solute charge 
distribution p(r). We denote the one spherical boundary in the problem, which separates the protein and 
solvent, by &, and remind the reader that in this problem we are only treating a single boundary. Inside the 
protein, the total electrostatic potential ipi(r) again obeys the familiar local- response dielectric theory with 
dielectric constant e pr otein: 



El — — V (^proteinj 

Di(r) = eproteineoEi(r) 
V • Di(r) = p(r). 

We denote the Coulomb potential due to the fixed protein charges as 



¥>mol 



Q 

E 

k=l 



tprotem | 



(12) 
(13) 
(14) 



(15) 



and the reaction potential due to the difference between the protein and solvent dielectric properties by 
tpre&c, the total electrostatic potential is 



^i(r) = ^ m oi(r) + v? IC ac(r). 



(16) 



In this nonlocal problem, we have a pure water solvent (no mobile ions) in which the displacement and electric 
fields are related nonlocally by a convolution with a dielectric function of the form £(r, r') = e(|r — r'Q so 
that 



D„(r)=e / £(r,r')E„(r')rf 3 r 
V-D ir (r) = 0, 



(17) 
(18) 
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and e(|r — r'Q is the Lorentz nonlocal function 



f (r, r') = e^r - r') + ^ ^ exp(-|r - r-|/A) ^ (ig) 

where e W ater is the bulk solvent dielectric constant (80 in the present work), is the short-range dielectric 
constant, here taken to be the optical dielectric constant 1.8, and A is an effective parameter that reflects 
the length scale associated with correlations between solvent molecules. At the solute-solvent interface 6, 
the usual Maxwell boundary conditions Eqs. |4] and [5] apply. By Eqs. [T7land[T8l the potential in the solvent 
must obey not the familiar Laplace equation but instead the hitegrodifferential equation 



V- f £(r,r')V^ n (r')dV 
Ju 



0, (20) 



the solution of which requires substantial calculation even for simple cases such as a sphere with central 
charge [TSJ [THJ HH1 HH1 HI] or a charge approaching a planar half-space }2"9l Hoi IT71 |4"5] . 

Hildebrandt et al. recently reformulated this nonlocal model as a system of coupled but purely local 
partial differential equations (PDEs) [55]. Similar simplification strategies have been demonstrated for 
modeling dispersive electromagnetic media |43| and plasticity |22| . Essentially, for a nonlocal relationship 
that takes the form of a Green's function for a known PDE, one may be able to introduce a new unknown 
potential whose gradient is the vector field resulting from the convolution (here Dn). Enforcing the original 
conservation law (here, V • D = 0) leads to an additional Laplace equation and then the original unknown 
interest and the additional unknown are coupled. For the Lorentzian model, the nonlocality resides in the 
second term of Eq. [TH] which is merely the Green's function of the Yukawa equation V 2 u(r) = A 2 u(r). Here, 
by introducing the auxiliary displacement potential ipu, one may write the coupled PDE system as 

VVW="P(r) (21) 
V^u(r) =0 (22) 

V 2 - j^j Mr) = - ^'ii(r) (23) 



with A = Ay eoo/es- The exact displacement boundary condition (Eq. [5]) is nonlocal and slow to compute, 
and so Hildebrandt [29] proposed the approximate boundary conditions 

^(r-) =^„(r+) (24) 
d _ d 

eoe P rotein^-Vi(r^) =^-^n(r^) (25) 
on an 

d d 

— Vn(r+) =e e M ^n(r h + ). (26) 

Different choices for boundary conditions are analyzed in more detail elsewhere, with model calculations 
suggesting that the impact on many calculations should be small compared to the overall differences between 
local and nonlocal models [5B] . 

For numerical scaling, it is useful to change variables by introducing the substitution 

* = — ( — V'H - eproteinVW J , (27) 

as discussed extensively elsewhere [28) . Then, defining 

, /I t^y i e protcin V DR \ ,„ f e protein TA y ^protein rrDR \ ^gmol / os \ 

= - I - - K A + K A I lfi mo l - — V A - V A I , (28) 

\ * ^solvent / \ too ^solvent / 
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the complete BIE system is 



A 

K L 



-vl 



tprotoin y 



D7? 



A A 



^protein "^A-f/ 





VII 




' b ' 













dn 














(29) 



where Vf R 



V L and similarly K? R = KY- K l . 



We omit the lengthy derivation and refer interested 
readers to Hildebrandt [28] . 

A point of great importance for fast numerical solution of Eq. [29] is that each non-zero block is a linear 
combination of the same boundary integral operators as are needed to solve Eq. 1111 As a result, the same 
fast BEM solvers used for local electrostatics in the LPBE model (e.g., fast multipole methods [41], pre- 
corrected FFT [3S], and the FFTSVD algorithm [H[5]) can be adapted easily to solve nonlocal electrostatics 
models [11] , Fast solvers allow the discretized linear system, which is dense in the sense that the number 
of non-zero entries grows quadratically with the number of unknowns, to be solved in linear or near-linear 
time. 



2.3 Eigenfunction Expansions of Boundary-Integral Operators on Spheres 

All of the boundary-integral operators of Eqs. [Til and l29l arc diagonalized by the surface spherical harmon- 
ics [3T]. Consequently, the boundary integrals of the form J F(r,r')u(r')dA' can be re-written as 

/oo +n . 
F{v,v')u{v')dA' F m(M)A£ m / Y^{9'^')u{9'^')dA' (30) 

n=0m=-n J 

where the (9, <f) are the angular coordinates for r, Y£(0, 4>) are the orthonormal surface harmonics, and \f lm 
is the eigenvalue for the n,m mode. Note that Eq. [30] represents a slight abuse of notation, in that the radii 
of the "source" and "destination" spheres are included only implicitly in the actual eigenvalues. Also, in 
this work, the eigenvalues of the relevant operators are independent of m, so we omit the second subscript 
in the remainder of the text. 

For a sphere of radius i?, the eigenvalues of the four "self-to-self" operators V L , K L , V Y , and K Y are 

^ = 2^1 ^ 

A " L = ~2(2n+l) (32) 
X yY = i(i K )R 2 j n (iK,R)h^(iKR) (33) 

\* Y = i(i K ) 2 R 2 /2 (j n (iKR)hW(iK,Rf) (34) 

where i = \/— T, j n {x) and hn {%) denote the spherical Bessel function and spherical Hankel function of the 
first kind, respectively, and the prime notation in Eq. l34l dcnotcs differentiation with respect to the argument. 

The Kirkwood problem also involves four Laplace boundary-integral operators that map between con- 
centric spheres. We demonstrate in the Appendix that the eigenvalues of these operators are 

A„ a,b = I ' /i,\n+i i _ (36) 




Xn b ' a = { , ^ ( a \n _i „ ^ n ( 38 ) 

2(2n+l) ' 
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3 The Boundary-Integral-Equation + Eigenfunction Approach 



3.1 Application to the Kirkwood Problem 

To simplify the coupled boundary-integral equations, we introduce the spherical-harmonic projection opera- 
tor Y * , which maps a function defined on a sphere (i.e. in angular coordinates) into the expansion coefficients 
in the basis of surface spherical harmonics, which is complete and orthonormal. Similarly, the operator Y 
maps a vector of expansion coefficients in the basis of surface harmonics to a function on the sphere. 
The non-zero blocks of the matrix in Eq. [TT] can be simultaneously diagonalized as 
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(i)lH=o, Dy> = A^ w | fi=a , = A£ w |i{ =0 , and = A; (i) |R= a , where n(i) denotes the degree 

associated with the ith eigenmode. Expanded in the surface harmonics, the unknowns of Eq. QT]are written 
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and projecting the right-hand side similarly, we obtain the surface-harmonic analogue to Kirkwood's result: 
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Note that this representation does noi diagonalize the entire operator, but does decompose the reaction 
potential in the protein into the individual harmonics. 

An algorithm to solve the Kirkwood problem using the BIE/cigenfunction approach is therefore structured 
as follows. For each mode i to be solved (up to a desired order), one first computes the projection of the 
solute charge distribution onto the ith solid spherical harmonic (i.e. one computes the appropriate multipole 
expansion coefficient). Then one calculates the ith eigenvalues for the boundary integral operators to set up 
a linear system of equations with four unknowns, and solves for the ith expansion coefficient of the reaction 
potential. The reaction potentials at all desired locations is then easily computed. 



3.2 Application to Nonlocal Electrostatics 

We now derive our main result — the exact analytical solution of nonlocal electrostatics for a spherical solute. 
The 3-by-3 block operator of Eq. [29] can be decomposed as 
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where again Y* projects from a distribution on the sphere surface into an expansion in surface spherical 
harmonics, Y represents the harmonics themselves, and the matrices D^ 1 are all diagonal. The entries of 
the matrices are simply the appropriate scaled sum of the operator eigenvalues: e.g., D 
where n represents the degree associated with the ith eigenmode. Projecting both sides of Eq. 

D« D(») 

£)(6) £>(7) 

The ith entry of the projected form of Eq. [28] is therefore 





fin 




' b ' 




d(f>\\ 























- 1 \Ka 

— 2 A n(i)' 

one obtains 



(43) 



X' 



(i) 



Cprotein *Kj? R 

~ A n(i) 
fc solvcnt 



Pn 



: protcin >Vj 



n(i) 



^protein yV, 
^solvent 



DR 

00 



dip 



mol 



On 



(44) 



Again, solving analytically for each coefficient (p n u-\ independently provides the desired expansion (in surface 
harmonics) of the potential at the protein-water boundary. These coefficients are readily converted to 
the solid harmonics to obtain the potential inside the sphere. The analytical nonlocal model has been 
implemented in MATL AB and is available as Supplemental Information [TJ [14] . 

It may be verified that in the limits A — > and A — > oo, the analytical solution converges to the appropriate 
local-response models; sec Figure [5] for the example of a sphere with a single central charge, which is known 
as the Born ion. As a more challenging validation, we have used the nlFFTSVD fast BEM solver [TT] to 
compute the solvation free energy of a single +le charge situated at (0, 0, 6 A) inside a sphere of radius 
8 A centered at the origin, and plotted the convergence of these results to the solvation free energy computed 
analytically (Figure [3]) . This test case is challenging because it lacks the spherical symmetry of the Born- ion 
test case, and in fact BEM simulations require finer discretization for charges close to the surface [33] . 

The required spherical Bessel and Hankcl functions have been computed using the algorithm proposed by 
Cai |19] , and their derivatives were calculated using well-known recurrence relations [53] . Using numerically 
stable implementations of the Bessel functions and their derivatives is of utmost importance. For large sphere 
radii and charges approaching the surface, large cancellations in naive implementations of the projections 
causes the observed value of the solvation energy to diverge as the order of the calculation is increased. 
Additionally, very large and very small values of A arc problematic for the calculation of the Yukawa-operator 
eigenvalues, and suggest that further research into their accurate computation is warranted. 



4 Results 

Hildebrandt et al. have suggested that A ~ 15 — 24 A provided an excellent fit to experimental data for 
monatomic cations [29] [28]. However, other factors, especially nonlinearities such as dielectric saturation [25], 
may play important roles in ion solvation and charge burial, it is important to understand the dependence 
on A. Figure 2] contains plots of nonlocal-model electrostatic solvation free energies for monovalent cations 
of varying radii. The nonlocal-model free energies are clearly sensitive to A in the range 1 to 10 A, but 
much less so outside that range. Therefore, although ion solvation free energies therefore provide a clear and 
intuitive demonstration of the impact of nonlocal response, the insensitivity outside of A > 10 A suggests that 
that such data should be used with caution when performing more detailed parameterization; in particular, 
we emphasize that it is impossible to fit the ion radii as well as A simultaneously because the problem is 
underdetermined. More extensive calculations are needed to calibrate the new model, and are underway. 

Analytical solutions for simple geometries also allow fast determination of the reaction potential through- 
out the whole system. More thorough visualizations of solvent response may offer new insights into the 
empirical, seemingly application-specific definitions of the protein dielectric constant [551 HH] : including for 
example why values of e pro t e in much larger than experimental estimates |24] are often needed to obtain accu- 
rate calculations of pK a shifts in proteins [20] ■ To illustrate the fundamental differences between local and 
nonlocal theory, as well as the computational advantage of having a fast analytical model for visualization, 
we plot the reaction potentials for both simple and complicated charge distributions as we vary key model 



9 



2 

10 □□□□□□□ nn Dnnn . nOOOOOOOOOOOOOOOOOOOOOOOO 

f ° a "u e o° ° 

3 o 
r o° 

2 m o n Q 

< 10 ■ ° a 

o 

8 i ° □ 

1Q- 1 T 



o A. = limit 



□ 
□ 

□ 

□ 

□ 



10 



□ A, = oo|imit ° 

"2 , I , 1 , □ 



-2 2 4 

10 10 10 10 



X (Angstrom) 



Figure 2: The analytically computed solvation free energy for a sphere with central charge (Born ion) 
converges to the correct local-response limits as the nonlocal length-scale parameter A approaches or oo. 

parameters: the protein dielectric constant in the local theory, and the effective length scale A in the nonlocal 
model. 

Figure [5] contains plots of the reaction potential induced by a single +le charge in a protein-sized sphere 
of radius 24 A, where the charge is situated 2 A from the dielectric boundary. The reaction potential for 
local-response models is shown in (a) and (b), with e pro tem = 2 in (a) and e pro tein = 4 in (b). Nonlocal- model 
results are plotted in (c) and (d); for both nonlocal calculations, e pro toin = 2, with A = 1 A in (c) and 
A = 10 A in (d). For comparison, all potentials are plotted according to the same color scale. Adjusting 
^protein from 2 to 4 in the local model leads to a qualitative global shift in the reaction potential. On the 
other hand, nonlocal response presents relatively small overall changes, even though A varies substantially. 
For a single +le charge buried deep within the protein at (0, 0, 10 A), the reaction potential is smaller in 
magnitude, which means that the qualitative shift for increased e pro tein can be seen more easily (Figure [6]). 

These qualitative differences are meaningful for the types of complicated charge distributions found in 
proteins as well. To illustrate this point, we use as an example the protein bovine pancreatic trypsin inhibitor 
(BPTI). We model the charge distribution by taking the atomic coordinates from the Protein Data Bank 
(accession code 3BTM [57]) and assigning atomic charges using the PARSE [53] force field. Figure [7J contains 
plots of the resulting reaction potentials; the results for each subfigure are computed using the same model 
and parameters as used for the corresponding subfigure of Figure [3J Together, these results suggest that 
future nonlocal studies should investigate charge-charge interactions in more detail, especially contrasting 
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Figure 3: Relative error for numerical simulations of the nonlocal model using BEM, as a function of the 
number of unknowns in the discretized problem, for a 8-A-radius sphere with a single +le charge situated 
2 A from the sphere surface. 

the fields induced by buried and surface-exposed charges. 

We would like to emphasize the substantial difference in speed between numerical and analytical methods. 
For the BPTI test problem, a low-resolution numerical simulation using the highly optimized, linear-scaling 
boundary-element method (BEM) code nlFFTSVD [TT| — one of the fastest numerical implementations of the 
nonlocal model — requires approximately 12 minutes on a 2012 MacBook Air. The unoptimized MATLAB 
implementation of the analytical approach, in contrast, is more than 45 times faster, requiring less than 17 
seconds. 

5 Discussion 

The shortcomings of local electrostatics continue to motivate new models, but often the practical compli- 
cations of numerical simulation slow their testing and improvement. To accelerate studies of the promising 
Lorentz nonlocal model 1301 HH) j we have derived the exact analytical solution for a spherical solute 
containing an arbitrary charge distribution. Our approach uses Hildebrandt's boundary- integral equation 
(BIE) formulation |30| and the analytically known eigendecompositions of the associated boundary-integral 
operators. Calculations demonstrate the method's correctness and that solvent screening of charge-charge 
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Figure 4: Dependence of the electrostatic solvation free energy for a +le point charge centrally located in a 
sphere, as a function of the sphere radius and the nonlocal parameter A. Here e pro tein = 1, ^water = 80, and 



interactions are markedly different in nonlocal and local theories, even when the protein dielectric constant is 
adjusted. Fast analytical models enable rapid visualization of electrostatic fields, and thus facilitates efficient 
exploration of the new model's implications and qualitative differences from existing theories. 

The BIE-cigenfunction strategy represents a novel alternative to matching potential expansions and may 
be useful in other areas of mathematical physics. To illustrate the method's generality, we have also derived 
the solution to the Kirkwood two-boundary problem for local electrostatics, which has furnished many 
insightful physical studies and model approximations even though proteins clearly take shapes much more 
complex than spheres. It should also be noted that a similar algorithmic approach can simplify calculations 
involving matched solid-harmonic expansions; that is, instead of tedious, error-prone algebraic manipulations 
to obtain the desired expansion coefficients, one could set up the small linear systems for each mode and 
allow the computer to do the arithmetic. 

The present work enables studies of the nonlocal model to be conducted rapidly for simple model systems, 
obviating the need for more complicated and slower numerical calculations [30l [57j [TTJ HH] . To encourage 
further tests of nonlocal models, the Supplemental Information [TJ [2] includes a MATLAB implementation 
of the analytical approach. As described in earlier work on nonlocal electrostatics, boundary conditions 
represent a subtle issue that warrants detailed study [2H1 \E7\ [23] , and fast calculations on spheres will allow 
a simple way to test improvements. Our results also provide a useful way to test numerical simulations 
of nonlocal electrostatics on nontrivial systems, e.g. models of finite-sized solutes with complicated charge 
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Figure 5: Reaction potential (in kcal/mol/e) induced in a sphere of radius 24 A by a single +le point 
charge situated 2 A from the boundary, for different local and nonlocal models. All potentials are plotted 
on the same color scale, and for all models, e W atcr = 80. (a) Local-response model with e pro tein = 2; (b) 



local-response model with e pro tcin = 4; (c) nonlocal-response model with e prc 
nonlocal-response model with e pro t C m = 2 and A = 10 A. 



2 and A = 1 A; (d) 



distributions. 

Future work will address the development of fast analytical approximations similar to recent Generalized- 
Born (GB) models [51 j or BIE approximations [12| . Second-kind boundary-integral formulations may offer 
substantial advantages for such approximations [40] . and Fasel ct al. have recently presented a purely 
second-kind formulation of the nonlocal model |23j . An extension of our approach to the Fasel formulation 
is therefore of significant interest. One extension to the present work might be to account for the fact 
that many proteins can be reasonably well modeled using ellipsoids (see, for a recent example in electrostatic 
theory, [50] )■ It is possible that one could use a similar approach to derive an analytical solution for ellipsoidal 
geometries as well; the eigendecompositions of the Laplace boundary-integral operators for ellipsoids are 
known, for instance [2j [3[ I36| , though corresponding results for the Yukawa integral operators do not appear 
to have been published. We also note that even for the sphere, computing the eigenvalues of the Yukawa 
integral operators is numerically challenging, and should motivate the development of improved algorithms. 
Recent work on computing the ellipsoidal harmonics found similar challenges |13| . and the present work 
has uncovered a second compelling example of how molecular biophysics poses novel challenges for more 
fundamental research in applied mathematics and numerical analysis. 
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Figure 6: Reaction potential (in kcal/mol/e) induced in a sphere of radius 24 A by a single +le point charge 
buried 14 A from the boundary, for different local and nonlocal models. All model parameters are the same 
as the corresponding plots in Figure [5] 



Appendix: Eigenvalues of the Laplace boundary-integral operators 
for concentric spheres 

We first address the single- and double-layer operators that map from the inner sphere (radius b) to the outer 
(radius a). For the single-layer operator, let us expand a surface potential on the inner sphere in surface 
harmonics, i.e. 

i> Sb =Y^Sl m Ynifi,4>)- (45) 



and also expand the potential field in the region outside that sphere 



-{i»+l)yra 



(46) 



These two fields must agree on the surface r = b, and by orthogonality of the Y™ functions, we have 

V nm = S b nm b n+ \ (47) 
A similar surface expansion holds for the fields on the outer concentric sphere 

vx = E 5 «™w^ ( 48 ) 
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Figure 7: Reaction potential (in kcal/mol/e) induced in a sphere of radius 24 A by the charge distribution 
of the protein bovine pancreatic trypsin inhibitor (BPTI), for different local and nonlocal models. All model 
parameters are the same as the corresponding plots in Figured] See main text for full computational details. 



which may be matched to Eq. 05] to give 

V nm = a n+1 S% m . (49) 

Combining Eq. |47] with Eq. [49j we have 

] n+l 

<f» _ qb /m\ 

Because the eigenvalue for the single- layer Laplace surface operator on the inner surface is b/(2n + 1), we 
finally have that 

We derive the double-layer potential operators using an alternative approach based on Green's theorem. 
Consider again the expansion in spherical harmonics of the potential outside b from Eq. 1461 so that the radial 
component of the electric field is 

?± = Y,V nm {-{n+l))r-^Y™{0,ct>), (52) 



so the normal derivative of the potential at the inner surface b is 

dr 



^\ r=b = £ V nm (-(n + l))6-("+ 2 )y™(0, 0). (53) 
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Green's theorem allows us to write the potential at any point r with r = a > b as 

Hr) = + f d -^l^)dA' f G(r, i>)*tp.dA'. (54) 
Jb on J b dn 

Using Eq. 1301 and again relying on the orthogonality of the harmonics, we obtain 

a -(«+D = X<%-^) - \< b (-{n + l))^" +2 >; (55) 

v L 

substituting the known \ n a ' b from Eq. 1511 gives 

a-(" +1 ) = X^b-^ + ^±l a -("+D (56) 

" 2n + 1 v ' 

and finally 

*f-STi (:)"'■ (57) 

This result may be checked in the limit as a — > b, where Eq. [54] becomes 

^(r) = I^(r) + J ^p-^(r')dA' J G(r, r')^-dA'. (58) 



Analogous manipulations lead to the relation 



a -(n+l) = l a -(n+l) + \K b ~(n+l) + H±l cr (n+l) 

2 2n + 1 



(59) 



and thus we recover the self-surface result that = 2 (2n+i) • ^ ne e ig enva hies for the operators that map 
from the outer sphere to the inner one are obtained in very similar fashion using interior harmonics. For 
example, 

iP = Y,V nm r n Y™(6, ( j ) ), (60) 



and equating coefficients as before 
so that we have for the single-layer 



Vnm — ^ n S nm — ^ n S nm (61) 



The eigenvalues presented for these operators can be verified analytically using Green's theorem. 
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